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tween hydrophobic walls at partial hydration, with a wall-to-wall distance of about 0.5 nm. 
We analyze how the diffusion constant parallel to the walls, Dn, changes and correlates to the 
phase diagram of the system. We find a locus of Dm maxima and a locus of Dm minima along 
S^ isotherms, with lines of constant Dii resembling the melting line of bulk water. The two loci 

fa: 

of Dii extrema envelope the line of temperatures of density maxima at constant P. We show 
how these loci are related to the anomalous volume behavior due to the hydrogen bonds. At 
much lower T, confined water becomes subdiffusive, and we discuss how this behavior is a 
consequence of the increased correlations among water molecules when the hydrogen bond 
network develops. Within the subdiffusive region, although translations are largely hampered, 
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we observe that the hydrogen bond network can equilibrate and its rearrangement is respon- 
sible for the appearance of density minima along isobars. We clarify that the minima are not 
necessarily related to the saturation of the hydrogen bond network. 

Introduction 

Water displays many thermodynamic and dynamic anomalies. *~ 3 Although the origin of these 
anomalies can usually be traced back to the properties of hydrogen bonding between water molecules, 
their quantitative understanding is a goal to which a great deal of effort is being devoted. Experi- 
ments show that these anomalies are evident where liquid water is stable, but are stronger below the 
melting line where water is supercooled, i.e. metastable with respect to crystal ice. Bulk water can 
be kept in this supercooled state down to about 235 K at 1 atm. 4 The experimental limit of stability 
of supercooled water with respect to crystal ice defines the homogeneous crystallization tempera- 
ture Th(P), which reaches its lowest value of 181 K at about 2000 atm. In an attempt to rationalize 
these anomalies, several ideas have been proposed, including the stability limit (SL) conjecture, 5 
the liquid-liquid critical point (LLCP) scenario, 6 the singularity-free (SF) hypotheses,- and the 
critical point free (CPF) scenario. 8 All these stimulating ideas are consistent with the experimental 
properties of water, but hypothesize different behaviors below T^{P). Although these differences 
cannot be directly tested in experiments, their implications in the interpretation of water properties 
could be different, or could be relevant for other anomalous liquids.- - — It is therefore worthwhile 
to test these hypothesis in theoretical models. Many authors resort to computer simulations of de- 
tailed models of water (see for example 12 for a short updated list). However, this approach, both 
for molecular dynamics (MD)— or Monte Carlo (MC) simulations, 14 faces the problem of large 
computational times near T h (P), because water equilibration time increases in an exponential way 
for decreasing T. 

An alternative approach, which we follow here, is to consider a coarse-grained model of water 
that allows to perform, on the one hand, efficient MC simulations and, on the other hand, theo- 



retical calculations. In particular, we focus on a model for a water monolayer confined between 
hydrophobic walls.-^*^ The interest for this case comes from the fact that, under appropriate con- 
ditions, the formation of ice can be avoided in experiments with water under confinement. 17 ~ 19 We 
consider here only the case of confinement between infinite hydrophobic walls, while other kind of 
confinements, mimicking porous hydrophobic materials, have been considered in other works.— 

The coarse-grained model considered here allows to gain an insight into the properties of the 
scenarios that have been proposed for supercooled water. In particular, it is possible to show that 
within the framework of this model all the scenarios proposed for supercooled water differ only 
in the relative strength of the directional (covalent) and the many-body (cooperative) component 
of the hydrogen bond (HB). 21 When the many-body HB component is strong and the directional 
HB component is weak, the model recovers the CPF scenario, and shows that it coincides with the 
SL case. When the many-body HB component is zero, the model reproduces the SF scenario for 
any finite directional HB component. Finally, for intermediate values of the two HB components, 
the model recovers the LLCP scenario, with the liquid-liquid phase transition having a negative 
slope in the pressure-temperature (P — T) phase diagram. The LLCP occurs at positive or negative 
pressure depending of the relative strength of the many-body HB component with respect to the 
directional HB component.— Direct experimental evaluation of the relative strengths of the two 
HB components is not straightforward, but indirect evaluations are consistent with values that, for 
the model, would predict the LLCP scenario. 21 

Dynamic properties 

Among the many dynamics anomalies of water, we will focus here on the behavior of translational 
(self)diffusion constant D. In normal liquids D decreases when P increases at constant T, while 
water displays a large increase of D in a delimited region of P-T ,— "22 e. g., with an enhancement 
of about 60% at 243 K when P increases from 0. 1 MPa to 150 MPa. 23 The increase is observed up 
to about 200 MPa.- 2 ^ 

Computer simulations of bulk water with detailed models (see Ref.s— ~— and references therein) 



and lattice models (e. g.— ~— ) can reproduce, at least qualitatively, the anomalous behavior of D. 
By analyzing the microscopic structure of water molecules in the region of the anomalous increase 
of D, several authors have proposed a relation between the behavior of D and the structure of wa- 
ter. For example, Ref. 34 relates D to the configurational entropy, Ref. 35 associates the minima in D 
with a maximum in orientational order, and Ref. 36 shows, for a water- like isotropic potential, the 
connection of the anomaly in D, and in other quantities, with the density dependence of the entropy 
in excess over the entropy of the ideal gas. In particular, by classical molecular dynamic simula- 
tions it has been observed that the increase of P weakens the hydrogen bonds, and thus increases 
D. °' 29 This interpretation in terms of defects in the HB network can be extended to negative P.— 
A similar qualitative conclusion has been reached also by ab initio molecular dynamics showing 
that D is directly linked to network imperfections. 37 Nevertheless, a quantitative relation between 
the anomalous behavior of D and the microscopic structure of water is still missing. 

In confinement, experiments show controversial results. Without the aim of reviewing the rel- 
evant literature, we recall here, as examples, that the reduction of water D can be of up to two 
orders of magnitude between 260 K and 3 1 K, in hydrophilic NaX and NaA zeolites , 38 or that the 
viscosity of water between two hydrophilic surfaces, with < 1 nm interfacial separation, is seven 
orders of magnitude greater than that of bulk water at room temperature. 39 Other experiments dis- 
play that D decreases if the confinement increases, e. g. in MCM-41 with pore radius between lnm 
and 2nm— "^ or MCM-48 with pore radius of about lnm, 41 both slightly hydrophilic due to the 
oxygen atoms in the silica structure, with the first considered more hydrophobic than the second. 
Similar results were observed for channels of closed multiwalled hydrophobic carbon nanotubes 
with diameter between 2nm and 5nm. 44 Nevertheless, other experiments reveal an exceptionally 
fast mass transport for water confined in carbon nanotubes of about 2nm 45 and 7nm radius. 46 

Computer simulations have been performed to rationalize the different experiments. By MD 
it has been found a decrease of D with decrease of separation between two hydrophilic surfaces 
at nanoscopic distance for SPC/E water. 47 Moreover, it has been shown for the same model that 
the decrease of hydration largely decreases the mobility of water molecules near the surface of 



MCM-41 or Vycor, consistent with the interpretation that at low hydration the majority of them 
are bonded to the surface. 48 

In the case of hydrophobic confinement, the results are more controversial. MD of TIP5P 
water nanoconfined between hydrophobic smooth walls displays anomalous D,— but only in the 
direction parallel to the walls. 50 D results to be two orders of magnitude smaller than in bulk and 
with the anomaly occurring in confinement at lower T than in bulk.— A similar large decrease of 
mobility has been reported at ambient conditions for SPC/E water between two large hydrophobic 
graphite-like plates for separations below 1.3 nm. 51 Nevertheless, first-principle MD simulations 
of the same model in similar conditions show that the diffusion of water molecules become faster 
under confinement, possibly due to weaker HBs at the interface. 52 A similar controversy is reported 
for simulations of water in carbon nanotubes, with radii below 1 nm.- — 

Our approach 

Our approach to the problem is to study by MC the local dynamics of a coarse-grained model of 
confined water that will be defined in detail in the following. The results described in the previous 
section for this model have been derived by free-energy calculations within mean-field approach 
and efficient MC simulations. In particular, a cluster MC dynamics allows to easily equilibrate the 
model at any T and P. 56 This specific MC approach is based on a mapping of the thermodynamic 
model into a geometrical problem, using an appropriate extension of the correlated-percolation 
approach.— — Nevertheless, as an alternative it is possible to adopt a local MC algorithm with the 
aim of studying the dynamic behavior of the model and compare it with the experimental behavior 
of water M& 

This kind of study is useful for systems at equilibrium, although approaching a glassy dynam- 
l cs 59,63,64 anc | ma k es ^e plausible assumption that, at a given T and P, the MC time step can 
be converted into real time units by a factor that does not depend on time. This assumption is 
consistent with the mode-coupling theory (MCT), according to which the long-time relaxation (a- 
relaxation) dynamics should be independent of the microscopic dynamics, 65 as tested in Lennard- 



Jones mixtures. 66 However, the comparison must be performed with caution, because the time 
conversion factor for a given observable could depend on both T and P, when the interactions are 
non-isotropic. 67 This dependence can be estimated by comparing MC results with experiments, as 
done by Mazza et al. in Ref. 68 for the model under consideration here. The comparison shows 
a linear relation between the logarithms of time in real units and time in MC step at different T 
and constant i>, 68 corresponding to a power-law relation between the two times. This relation can 
be easily understood in the context of the present coarse-grained model, at least at low T. At low 
enough T, both the experimental time and the MC time follows a generalized Arrhenius relation 68 
where the characteristic r -dependent activation energy corresponds to the energy needed to break 
a HB. The power-law relation between the two times turns out to be a consequence of the choice of 
the model parameters that implies a HB energy lower than in experiments. Therefore, by choosing 
appropriate model parameters the two times would be directly proportional, with a proportionality 
factor that would not depend on T at constant P. 

With this caveat in mind, we perform here an extensive MC study of the thermodynamic and 
dynamic behavior of the coarse-grained model, presented in Ref., 15,16 for a water monolayer con- 
fined between hydrophobic walls. In particular, we consider the case in which water molecules 
can diffuse. We identify the region where diffusion has an anomaly behavior, finding maxima and 
minima of the diffusion coefficient at fixed T, and we discuss how this anomaly is related to other 
anomalies of water. 

In particular, thanks to the feature that the model can be tuned from one scenario to another, 21 
we can test if the diffusion anomaly is related to some specific scenario for the thermodynamics of 
the system at much lower T. Specifically, we consider the LLCP and the SF scenarios. We do not 
observe any relevant differences in the region where the system displays diffusion anomaly. 

We finally investigate the very low T dynamics, observing a subdiffusive regime and the ap- 
proach to a glassy state at any P. Under these conditions, we observe density minima, resembling 
recent experiments 69 and MD simulations. 70171 In the present model, the density minima appear 
in the vicinity of the glassy state and as a consequence of the breaking of HBs to rearrange water 



molecules for a better matching with the local order. In particular, we observe that if the model 
does not take into account the many-body interaction, the density minima are largely reduced. 

Coarse-grained model of a water monolayer 

We consider a water monolayer confined between two hydrophobic smooth walls with double 
periodic boundary conditions. The walls are separated by a distance of about h = 0.5 nm. The 
system is considered at constant number N of water molecules, constant T and constant P, leaving 
the volume V free to change. For a TIP5P-water monolayer confined between hydrophobic walls 
it has been observed that, depending on the separation h, water is liquid or forms a quasi-two- 
dimensional square ice for temperatures ranging between 260 K and 300 K, and negative lateral 
pressure. 49 - 72 The square symmetry is a consequence of the distortion imposed by the walls to the 
tetrahedral HB network that would otherwise form in bulk water. Consistently with these findings, 
we divide the available volume V into jV square cells, each with volume v = V / jV , square section 



of size yjv/h and height h, and hydrate the system with N < ^V water molecules. 

In our coarse-graining we assume that each cell can at most host one water molecule. There- 
fore, if N = jY , then each cell has one molecule and the system is homogeneous in density. 
If N < jY , some cells are empty. To each cell we associate an occupation variable n ; = 0, 1 
(i — 1,2,... jV) if it is vacant or occupied, respectively. Having only one water molecule per cell 
between the walls, for the sake of simplicity, we reduce the description of the monolayer to a 
two-dimensional system. 

The water-water interaction is decomposed into three terms. The first accounts for all the 
isotropic contributions, including short-range electronic orbitals repulsion and van der Waals at- 
traction, and is represented by a Lennard- Jones potential truncated at a hard-core diameter r$ 



for r <tq 
for r > ro, 
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where e is the interaction energy and r is the distance between two molecules. In the coarse-grained 
representation, r is the distance between the center of occupied cells. The hard-core diameter 
tq is introduced to simplify the implementation of the model and our tests show that the results 
do not depend on its existence. We set tq = 2.9 A, the water van der Waals diameter, — ^ and 
e = 5.8 kJ/mol, consistent with the value 5.5 kJ/mol of the estimate of the van der Waals attraction 
based on isoelectronic molecules at optimal separation.— 

The second term of water- water interaction accounts for the directional (covalent) component 
of the HB formation, 1 ^LL with a characteristic energy that we set J = 2.9 kJ/mol. To account for the 
directionality we adopt a geometrical definition in which the HB breaks if OOH > 30°. Therefore, 
only 1/6 of the orientation range [0, 360°] in the OH-0 plane is associated with a bonded state. We 
therefore associate to each molecule i a bonding index cy G [1,2, . . .q] with a discrete number of 
states q describing the bonding state with a neighbor molecule j, and choose q = 6 to account for 
the entropy loss associated with the formation of a HB. Due to the square symmetry, each molecule 
has four neighbors and four bonding indices Oy. Therefore, each molecule has q 4 = 6 4 = 1296 
possible bonding states. To form a HB between two molecules in two occupied neighboring cells 
i and j (hence mtij = 1) we assume that the two facing bonding indices oy and Oji are in the same 
state, i.e. 8 ai . i0ji = 1, with 8 Chb = 1 if a = b, and 8 a ^ = otherwise. Therefore, the directional term 
of the HB can be expressed as 

Hj = -JN HB , (2) 

where 

N H B=Y, n i n jdou>°ji ( 3 ) 

is the number of HBs and the sum runs over all nearest-neighbor cells. 

The directional interaction of HBs leads to local reduction of density as a consequence of the 
reduced nearest neighbors with respect to close molecular packing. This property is at the origin 
of the macroscopic density maximum that occurs at 4°C at ambient pressure, a temperature below 
which the number of HBs per molecule is about 3. 78 The effect can be observed in experiments 
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mainly as a change in the structure of the second shell of water molecules.— 1 ^ - — Nevertheless, 
to include the effect in a tractable way in the coarse-grained model, we follow Ref . 7 and consider 
that each formed HB leads to a small increase of volume vhb per molecule, where vhb/vq = 0.5 is 
the average density increase from low density ice Ih to high density ices VI and VIII and vo = hr^ 
is our approximation to the van der Waals volume of a molecule. Therefore, the total volume 
occupied by water is 

V w =Nv + N hb vhb- (4) 

Note that the increase vhb corresponds to a larger volume per molecule, but not to a larger separa- 
tion r between molecules, hence it does not affect the radial term in ??. 

The last term we include in the water interaction energy is a many-body (cooperative) interac- 
tion among HBs, that favors specific values of the probability distribution of 0-0-0 angles— 1 ^- 
(see also Ref. 25 for a brief description of the quantum origin of the cooperative interaction). Fur- 
thermore, the probability distribution of 0-0-0 changes when comparing bulk and confined wa- 
ter, 82 showing the disappearing of the fifth interstitial neighbor in the confined case, and a shift of 
the maximum of the distribution toward 90° at low T, consistent with the symmetry chosen for the 
coarse-graining in our model. To account for this cooperative interaction, we include in the model 
the term 

#Coop = —Ja2^ n i H 8 0ik ,c iV (5) 

' (*,/)* 

where (k, /),■ stands for the six different pairs of four bonding indices of a molecule i, and J a is 
related to the energy gain when the bonding indices order in the same state, with the maximum 
gain — 6J a per molecule, corresponding to the fully ordered case. 

As shown by Stokely et al.,— by setting the parameters £, J and vhb to finite values and 
tuning the parameter J , it is possible to reproduce all the scenarios that have been proposed for 
supercooled water. Here we consider two cases: (i) J a = 0.29 kJ/mol, corresponding to the LLCP 
scenario, and (ii) J a = 0, corresponding to the SF scenario. 

The choice in (i) is consistent with the experimental measure of HBs in ice Ih, approximately 



3 kJ/mol stronger than in liquid water. 83 If we entirely attribute this increase to the cooperative 
component of the HB, 84 we find J a ~ 0.25 kJ/mol for the two molecules forming the HB (each 
with an energy gain in absolute value 6J a ), which is consistent with our choice. 

The water-wall interaction is represented by a hard-core exclusion. Although the interaction 
of water with a hydrophobic, infinitely large, object could have a small attractive van der Waals 
component or a soft repulsion, for the sake of simplicity we assume here that the main effect of 
the confining walls is to inhibit the formation of ice, as observed by Zangi and Mark for h < 0.51 

72 

nm.— 

Therefore, the enthalpy of the system at pressure P is 

H = U(r)+Hj+ //coop + PV W . (6) 

For a given occupancy ratio N/,yV, the state of the system is fully specified by the average number 
density N/V and the set of Oy. 

The Monte Carlo method 

We perform MC simulations in the NPT ensemble for a system partitioned into JV = 2500 square 
cells and an occupancy ratio N '/ ' JV = 0.75%, corresponding to N = 1875 water molecules. Since 
we allow for changes of the volume in the direction parallel to the walls, the control parameter P 
represents the pressure parallel to the walls. To test for size effects, we consider also JV = 400 
and jV = 1600 at the same 75% occupancy ratio, with TV = 300 and TV = 1200 respectively. Our 
results do not show any appreciable size effect among these three cases. Likewise, we do not find 
any significant differences for occupancy ratios between 75% and 90%. 

MC and mean-field results for this model at full occupancy— &1M&2.-22. were reported in pre- 
vious works i i£a22a2LSL£2j£2i22r-& An analysis of the differences between partial and full occupancy 
will be presented elsewhere. 

To generate equilibrium configurations we pick at random a cell i E { 1 , . . . , N} and an integer 
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j G {0, 1,2,3,4}. If j = 0, we choose at random one of the four neighboring cells of the cell i 
and, if empty, we displace into it the molecule i with probability given by the following ?? of a 
Metropolis algorithm. 

If j > 1 we set the bonding index a, 7 to any of the q possible states, independent of its original 
state, and accept the change with probability in ??. One MC step consists of 5N of these trials 
followed by a volume-change attempt, in which we select a new volume at random in the interval 
[V — 5V, V + 8V] with 8V = 0.5vo, and accept the move with probability in ??. Since we change 
the volume V in a continuous way, the volume per cell v and the distance between the center of cells 
r also change in a continuous way, as in an off-lattice model, despite the fact that the model has 
a fixed maximum number of nearest neighbors, equal to four, due to the square lattice symmetry 
adopted in the coarse-graining. 

From the new and old configurations we calculate AH = // new — H old from ??, and AS = 
— /Vfcgln(V new /V old ), where kg is the Boltzmann constant. We accept the new configuration with 
probability 

min{l,exp[-/3(AH-TAS)]}, (7) 

where /3 = \/(k B T). 

Following Ref. 68 we adopt real units to represent our results. The transformations to real units 
are based on rescaling the MC T, P, and t from experimental data for a monolayer of water ad- 
sorbed on lysozyme powder 68 and adjusting the MC p to experimental values around the maximum 
density at ambient pressure in a self-consistent way. These transformations are meant to give the 
order of magnitude of the calculated quantities. 

In the calculation of the Lennard- Jones interaction energy in ?? we test that there is no appre- 
ciable difference if we introduce a cut-off at the 9th neighbor for the maximum interaction range. 
To allow for a better equilibration of the system, we follow an annealing protocol along isobars, 
starting at high T, with jY cells randomly occupied by TV water molecules, each molecule having 
a random configuration of bonding indices, and with the total volume V = ^Vvq. We equilibrate 
each state point for 0.2 ms and produce data for 15 ms. 
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We calculate the coefficient Dm of water diffusion parallel to the plates from the Einstein rela- 
tion in two dimensions 

where r,(?) denotes the projection onto the plates of the position of molecule i at time t, and (■) 
stands for the average over all molecules i and over different values of to. 

To avoid correlations in the calculations of Dm and all the other quantities, we perform averages 
over blocks of To = 0.8 /is, sampled every 80 ps. We check that To approximately corresponds to 
the time needed for a molecule to reach its image points for most of the temperatures investigated. 
However, for the lowest temperatures, we use up to IOto as block-length for averages. 

Results and discussion 

Phase diagram 

At high T we find the gas phase, separated from the liquid phase by a first-order phase transi- 



tion ending in a critical point ( [figure][l][] 1 ). The liquid-gas critical point for the hydrophobically 
confined monolayer occurs at a pressure and temperature that are higher than that of bulk wa- 
ter, qualitatively consistent with the results of MD simulations for TIP4P water in hydrophobic 
confinement. 96 

By annealing the system from high temperatures, we find a discontinuous change in density 
along isobars ([figure] [2] []2). This change occurs at the spinodal temperature T^ G (P), that marks 



the limit of stability of the gas phase with respect to the liquid phase. The change is very large 
at low P and vanishes as the liquid-gas critical point is approached. At the critical point, by 
definition, the isobar has infinite negative slope. At P higher than the critical pressure, the min- 
imum slope of p(T) become finite and decreases in absolute value for increasing P. The locus 
of these minima corresponds to the locus of maxima of the isobaric thermal expansion coefficient 
a™ ax (P) = minx{— <9ln(pvo)/(9T|p}, with the maxima decreasing in value by increasing P. 
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Figure 1 : P — T phase diagram of a water monolayer nanoconfined between hydrophobic plates 
in which we emphasize how the diffusion constant Dm (z-axis and color scale) changes in relation 
with the thermodynamic behavior. The liquid-gas phase transition line (open diamonds) ends in 
the liquid-gas critical point (diamond symbol at highest T). The loci of isothermal maxima of Dm, 
D|V iax (solid squares), and minima, Dj" m (open circles), include the TMD line (solid triangles). 
Loci at constant Dm (e.g. the dashed line marked as Tso-D') resemble in their reentrant behavior 
the water melting line. Dm values above 0.12 A/ps 2 (gas phase) and below 0.03 A/ps 2 are not 
color-graded. State points with D|| <0.03 A/ps 2 are below the onset T (P) (dotted line) of the 
subdiffusive regime and their diffusion coefficient cannot be defined. For T < T (P), we observe 
a locus of density minima (asterisks), a locus of discontinuity in density (solid circles), a locus 
of maxima of Cp (dot-dashed line) as in a liquid-liquid phase transition ending in a critical point. 
Here we use e = 5.8 kJ/mol, J = 2.9 kJ/mol, J a = 0.29 kJ/mol, vhb/vq = 0.5 as parameters for the 
coarse-grained water model. 

By changing the simulation protocol across T^ G (P), i.e. heating the system instead of anneal- 
ing it, we find the typical hysteresis associated to a first-order phase transition, with the hystere- 
sis vanishing as the liquid-gas critical point is approached. All these results are consistent with 
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Figure 2: Density p as a function of temperature T along isobars at 0.02 GPa < P < 0.22 GPa 
separated by 0.02 GPa increments, as calculated by MC annealing. (Left panel) At low P, the 
discontinuity in p marks the gas-to-liquid spinodal temperature T^ G (P), vanishing at higher P into 
the liquid-gas critical point. The dashed line marks the temperature Tmd(P) of maximum density. 
For clarity we show only a selection of the simulated state points. (Right panel) Enlarged view 
at low T along isobars at 0.02 GPa < P < 0.30 GPa separated by 0.02 GPa increments. Model 



parameters as in [figure] [!][]! 
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previous findings for this model.— i£^i2L2£i£2r.£2i££i2£ j& Here we add the information about the 
diffusion constant. 

Diffusion maxima and minima 

Upon crossing the liquid-gas phase transition, the diffusion constant D» displays a discontinuous 
change that vanishes as the critical point is approached. Above the critical point, Dm changes in a 
continuous way, as expected in a one-phase region. 
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Figure 3: Diffusion coefficient Dm as a function of pressure along isotherms for (from bottom to 
top) T =483 K, 588 K, 658 K, 693 K, 762 K, 832 K and 1006 K. For T > 832 K, D\\ is monotonic 
with P, while it has maxima D|| iAX and minima Djj 11 " for 588 K< T < 832 K, and a maximum 
for T =483 K. The pressure p D ' MAX for DJ| 1AX (dot-dashed line) decreases for increasing T and 
converges to the pressure p D ^ min for Dj? 1 " (dashed line). p°- mm decreases for decreasing T and, 



possibly, becomes negative. Model parameters as in [figure] [1][]1 
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In the liquid phase for T > 832 K, D\\ decreases monotonically for increasing P, but has 



an anomalous non-monotonic behavior at lower T ( [figure] [3] []3l. This behavior resembles the 
known diffusion anomaly of bulk water. In particular, we find that Du decreases below a maximum 
jjMAX f Qr d ecreas i n g pressure at P < 0.2 GPa, as observed in bulk water. 24,25 The pressure at 
which we find D|V IAX , p D ~ MAX ^ increases with decreasing T. At lower P, D» reaches a minimum 
Djj 1111 at a pressure p D mm that decreases with decreasing T and eventually becomes negative for 



r<588K([figure][l][]l). 



As a consequence of the occurrence of D^* and Dlf™, the state points with the same value of 
Dii forms lines in the P-T phase diagram that are not monotonic as a function of P. Therefore, these 



lines at constant D», or iso-Dy lines ( [figure] [1][]1 ), have a positive slope in P-T phase diagram 



for P > P D - MAX and for P < P D ~ min , but a negative slope for pD-MAX > p > p D -mm It is 
interesting to observe that this change of slopes in P-T plane resembles the change of slope of the 
bulk water melting line, at least for P > p D_imn . Although the present model does not include any 
representation for the crystal, because we make the hypothesis that crystallization is avoided under 
the conditions considered here, our finding of the change of slope of the iso-Dii lines suggests that 
the shape of the melting line in nanoconfinement could resemble that of bulk water. Moreover, our 
observation hints that its shape would be mainly determined by the slowing down of the dynamics. 
Finally we observe that the temperatures of D ] u lAX and Djj 1111 are higher than those found in 
confinement and bulk with MD simulations of TIP5P- water. 49 This difference could be consistent 
with the higher values observed here for the liquid-gas spinodal line. On the other hand, the 
pressure at which we observe the onset of Dy 1AX is consistent with experimental results for bulk 
water.— &■ 

Density maxima and expansivity minima 

By decreasing T along isobars, we find maxima in density at a temperature Tmd(-P) f° r P ^ 
0.2 GPa, as observed in bulk water— ^ ( [figure][2][]2[ ). Above P ~ 0.2 GPa the density increases 



regularly for decreasing T, while at lower P and moderate temperature T < 7md(P), the density 
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of liquid water decreases for decreasing T . As a result of this anomalous density behavior the iso- 
baric thermal expansion coefficient (or expansivity) OCp(T) = — d\n(pvo)/dT\p becomes negative 
for T < Tmd(P), as observed in bulk water.- 

The temperature 7md(^) of maximum density has, in the P-T phase diagram, a shape that 
resembles the one observed in bulk and in nanoconfined water, 49 reaching a maximum temperature 
at about 50 MPa, similar to what is found in the TIP5P water model. 49 As expected for water, 35 
we find the Tmd locus within the region where the diffusion anomaly occurs, delimited between 



P u - MAA (r) andP u - min (r) ([figure][l][]P. 



For T < Tmd(P) we find that the slope of p(T) along isobars increases by decreasing T 
( |[figure][2][]2 ), implying a more pronounced negative value for <Xp, consistent with experiments 



for bulk 97 and supercooled confined water.— ^ The slope decreases at lower T, consistent with the 
occurrence of a minima in 0Cp, expected below 240 K from experiments in bulk water at ambient 



99 

pressure. y 



Subdiffusion 



Before describing in detail our findings about the density at lower T, we observe that for T < 275 K 
at low P = 0.02 GPa, within our simulation time ~ 0.01 s, the system does not reach the diffusive 
regime, i.e. the mean square displacement (|r,(/ + ?o) ~ r i(t)\ 2 ) is no longer proportional to the 
time t <\ [figure] [4] []4[). We find that the long-time behavior (t > 10 6 ps) is well described by the 



subdiffusive relation Ar 2 ~ t a , with a = 0.7. 

Subdiffusive behavior is observed in experiments for water hydrating mygloblin, at a hydra- 
tion level corresponding to approximately one water monolayer. 10 ° The experiments show water 
subdiffusion at 320 and 300 K, with an exponent a = 0.4. 10 ° This subdiffusive behavior (also 
called "anomalous diffusion") has been rationalized by several authors by means of simulations of 
water, both in inorganio 4 ^ 1 ^- and organic confinement,-^, - -^ 5 - with a exponents varying between 
0.96 103 and 0.45 ±0.05. 101 The proposed rationale is that the subdiffusive behavior is due to the 
heterogeneity of the surface and of the water-surface interaction. Nevertheless, this interpretation 
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Figure 4: Mean square displacement Ar 2 = (|r,-(f + fo) ~ r i(t)\ ) as function of time t. In log-log 
plot, the diffusive regime Ar 2 ~ t has a slope 1 (parallel to the line without symbols). We calculate 
the diffusion coefficient D\\ from ?? within the diffusive regime. (Left panel) At P = 0.02 GPa, for 
T < 275 K we find that the system does not reach the diffusive regime and the long-time behavior 
(t > 10 6 ps) is well described by the subdiffusive relation Ar 2 ~ t a , with a = 0.7. (Right Panel) 
At higher pressure P — 0.1 GPa the onset of subdiffusive regime occurs at lower T with respect to 
thei 5 



0.02 GPa case, for T < 205 K. Model parameters as in [figure] [1][]1 
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does not apply to our case, where the surface is by definition homogeneous and flat and the water- 
surface interaction is only due to volume exclusion. In our case the subdiffusive dynamics is, 
instead, originated by the increasing correlation among the water molecules that will be discussed 
in the following subsections. 

At higher pressure, P = 0. 1 GPa, we find that the onset T (P) of the subdiffusive regime occurs 



at about T < 205 K, i.e. at lower T with respect to P = 0.02 GPa ( [figure] [4] []4>. Therefore, 



within this range of P, the temperature T (P) is correlated with Ti so _£>(P) of iso-Di lines, being 



Tiso-D(^) — T (P) ~ constant. This is no longer true at P > 0.22 GPa ( |[figure][l][]l[ ) and can be 
understood in the framework of this model, because, as we will show in the next subsections, for 



P > 0.22 GPa the number of HBs vanishes at low T ( [figure] [6] [] 6 ). 



Density minima: relation with the cooperativity and the slow dynamics of the 
HB network 
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Figure 5: Density minima for the coarse-grained model of a water monolayer confined between 
hydrophobic walls. (Left panel) The locus of temperature T m [ n p (circles) of minimum density as a 
function of pressures P follows approximately a quadratic curve (continuous line) in the P-T plane, 
that extrapolates to about 201 K for atmospheric pressure. (Right panel) The density minima p m i n 
(circles) display an approximate square-root dependence on pressure P, whose extrapolation for 
atmospheric pressure is about 0.99 g/cm 3 . Model parameters as in| [figure] [!][]! 



Although for T < T (P) our MC simulations become subdiffusive, we find that we can equili- 
brate the HBs dynamics within our simulation time for temperatures as low as 190 K at 0.02 GPa, 
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Figure 6: Number of hydrogen bonds hhb to which a molecule participate, as a function of temper- 
ature T for pressures P from (top to bottom) 0.02 to 0.3 GPa in increments of 0.02 GPa. Inset: At 
low T for P from (top to bottom at 200 K) 0.04 to 0.12 GPa we find that hhb reaches a maximum 
value of ~ 3.75 at temperatures that coincide, within error bars, with those of the density minima 
p m i n at the same pressure. Model parameters as in| [figure] [!][]! 
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Figure 7: Comparison of the density behavior at low T for the model with parameters as in 
[figure] [l][]l] for P = 0.06 GPa (full squares) and P = 0.08 GPa (full circles), and with J a = 0, and 
the other parameters unchanged, for the same pressures (open squares and open circles, respec- 
tively). We find that for J a = the density minima is not detectable within our resolution. 
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Figure 8: As in [figure] [7] [] 7 but fori 3 = 0.02 GPa, with parameters as in [figure] [!][]! (squares) 



and with J a = (circles). The density minimum occurs also when J a = 0. 

or 163 K at 0.12 GPa. Specifically, we find that the relaxation time of the bonding indices C7y, 
related to the formation of the HBs, is of the order of 4 ns for these state points, while it exceeds 
our simulation time at lower T , e. g. at about 170 K for P = 0.02 GPa. 16 

In the region of state points where we can equilibrate the system, but close to the lowest well- 
equilibrated temperature, we observe a minimum in density along isobars with P < 0.12 GPa 



([figure] [2] []2). This result resembles the experimental density minimum for water confined in a 
nanoporous silica matrix MCM-41 with a pore diameter of 1.4 nm found by Mallamace et al. 69 

From our simulations for P > 0.02 GPa, at atmospheric pressure we extract a density mini- 
mum of about 0.99 g/cm 3 at about 201 K ([figure] [5] []5 ) not too far from the experimental value 
0.940 ± 0.003 g/cm 3 at about 203 ± 5 K and atmospheric external pressure. 69 Although our simple 
quadratic extrapolation predicts a value for p m [ n that is larger than the experimental, our data give 
an extrapolated T m [ n p at atmospheric pressure consistent with the results of the experiments 69 and 
comparable to those from MD simulations of TIP5P-E water. 70 

It must be noted, however, that the experimental results for confined water are controversial. 19 
Nevertheless, the controversy is mainly about the experimental measurement of the effect, and not 
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about the effect, because it has been observed that the existence of a density minimum in water is 
a necessary consequence of the existence of the low-P branch of the TMD line. 71 In particular, it 
has been proposed that the locus of p m i n corresponds to saturation, or the maximal ordering, of a 
network of water molecules with a random tetrahedral local arrangement. 71 

Our results, however, lead to a different explanation. We calculate the number of HBs hhb in 
which a molecule participates, defined as hhb = 2Nhb/N from ?? in such a way as to have four 
as maximum value for each molecule of the coarse-grained monolayer. First, we observe that our 
simulation results for the monolayer are consistent with experimental data for bulk water, with 
hhb — 0.45 at P = 0.25 GPa and T ~ 670 K, and with hhb — 2.2 at the lowest P at about ambient 
T, as reported in Ref. 81 (| [figure] [6] []6|). 



Next, we find that for P < 0.12 GPa the quantity reaches a maximum value of n#g ~ 3.75 at 
temperatures that coincide, within error bars, with the temperatures T m i n p of the density minima 
Pmin at the same pressure, and hhb decreases to ~ 3.5 below r mm p . Therefore, Nhb decreases and 
for ?? the density increases. 

To understand why uhb decreases below its maximum hhb — 3.75 at low T, we compare two 



cases corresponding to two different set of parameters of the model. The first as in [figure] [1][]1 
and the second, with J a = and the other parameters unchanged, corresponding to the LLCP and 
the SF scenario, respectively. 21 By comparing the low-T behavior of density for the two cases at 
intermediate P, we find that if J a = 0, i.e. the hydrogen bond is not cooperative, then the density 



minima is undetectable within our resolution ( [figure] [7] []7 ) 



However, at lower P we find that both sets of parameters give a detectable density minimum 



([figure][8][]8). Therefore, the cooperative term of the HB interaction is not essential for the 
occurrence of the density minima, but it emphasizes the minima at intermediate pressures. From 
this observation we understand that the explanation proposed by Poole at al.— can be applied to the 
case with J a = 0, in which the HB interaction does not include a cooperative (many body) term and 
the SF scenario is reproduced. In this case, the density increases for decreasing T at very low T, 
when all the possible HBs have been formed, generating regions of mismatching tetrahedral local 
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order. A decrease of T induces a small reduction of free volume per molecule, and a consequent 
small increase of density. 

Instead, when J a > 0, the cooperative interaction in ?? induces the breaking of HBs for de- 
creasing T, to allow the reorientation of a molecule and a better matching of local tetrahedral order 
at low T. As a consequence, the number hhb decreases from 3.75 to 3.5 inducing a large density 
increase. The high energy cost of this local rearrangement, i.e. the energy needed to break a HB, 
is at the origin of the high energy barrier for the process and the large increase of correlation time 
for the dynamics of the HB network in the vicinity of the locus of density minima. 

Relation of diffusion anomaly with different scenarios 

Several authors relate the diffusion anomaly in water to the presence of defects in the network of 
HBs. Here we show that the anomalous behavior of Dn in the coarse-grained model is not related 
to the many-body component of the HB interaction and, therefore, to the possible occurrence of 
the LLCP. The anomaly is, instead, due to the anticorrelation between volume and entropy, or to 
volume and energy, rooted to HB formation. 



We consider three different realizations of the coarse-grained model ([figure][9][]9). The first 



corresponds to the case presented in the previous section, with the parameters as in [flgure][l][]l 
For the second, we set J a = 0, leaving the other parameters unchanged. This case reproduces 
the SF scenario, 7 where density maxima occur and which has been shown to correspond to the 
vanishing-r limit of the LLCP. 21 Our simulations show that the occurrence of the anomaly of Dm is 



unaffected by this change of parameters ( [figure] [9] [] 9 ). Therefore, the absence of cooperativity in 



the HB dynamics is not relevant for the occurrence of both density maxima and diffusion anomaly. 
This can be understood for the clear separation between the temperature range at which the 
anomaly of Dn occurs and the temperature range at which water becomes subdiffusive. Only 
the latter regime corresponds to the temperature range for which the cooperativity has a strong 
influence on the dynamics, while it has no major dynamic effect at higher T. 



Next, we set vhb = an d leave the other parameters as in [figure][l][]l This case would 
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Figure 9: Effect of parameters vhb an d J a on the diffusion anomaly. We diffusion minima and 
maxima for both J a > and J a = 0. Instead, the non-monotonic behavior of Dn disappears for 
vhb = 0. 
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correspond to cooperatively bonding liquid with no density anomaly, i.e. with no anticorrelation 
between volume and entropy or volume and energy. In this case the change in the behavior of Dm 



is striking ( |[figure][9][]9[ ). The system has no diffusion anomaly, with Dm that decreases monoton- 
ically for increasing P as in normal liquids. 

Therefore, this result clarifies that the anomalous volume behavior due to the HB formation 
is directly related to the anomalous diffusion behavior. This conclusion, and the previous obser- 
vation that the many-body component of the HB interaction is not relevant in determining these 
anomalies, leads us to investigate how the volume available for diffusion, and the number of HBs 
are related to Dm, as discussed elsewhere.-^ 

Conclusions 

We study by Monte Carlo simulations a coarse-grained model for a water monolayer confined 
between hydrophobic walls. We consider a separation between walls about h — 0.5 nm that inhibits 
the formation of ice.— 

The model includes isotropic, directional (covalent) and many-body (cooperative) components 
of the HBs. Thanks to its coarse-graining, the model allows to study water in extreme conditions 
and, also, to check how each of the HB component affects its properties. Moreover, it makes 
possible to perform mean field calculations to compare with simulations results. 

We find gas and liquid phases, separated by a boundary line of first-order phase transitions 
ending in a critical point occurring at higher pressure and temperature, consistent with other models 
for hydrophobically confined water. 9 ^ We study the diffusion constant parallel to the walls Dm 
and find that it displays a line of maxima and a line of minima at constant T, as seen in similar 
confinement for other models. 49 Our analysis allows us to conclude that the anomalous Dm is a 
consequence of the anomalous volume behavior due to HB formation. In particular, the positive 
correlation between entropy and density, or energy and density, due to hydrogen bonding is the key 
element for the diffusion anomaly. It is worth reminding here that a similar result has been found 
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also for potentials with isotropic interactions and water-like anomalies 10 i 36 i 107 when the isothermal 
density dependence of the excess entropy, which is related to the total isothermal entropy by a linear 
function of the logarithm of the density, is considered. 

The difference with the present analysis is, nevertheless, threefold. First, for these isotropic 
potentials the HDL phase has less entropy than the LDL phase, implying a positive slope in the 
P-T plane for the liquid-liquid phase coexistence line as a consequence of the Clausius-Clapeyron 
equation. Instead, for water and the present model the HDL phase has more entropy than the LDL 
phase, hence the liquid-liquid phase transition has a negative slope in the P-T plane. 

Second, here we show that by setting the parameter that controls the increase of volume for 
HB formation, hence the positive correlation of density with entropy and energy, the anomalous 
diffusion behavior vanishes. Instead, for the isotropic potentials the vanishing of the anomalous 
diffusion behavior is controlled by the softness of the soft-core repulsion of the potential. 10111 A 
direct relation between these two results is not straightforward and could be interesting to investi- 
gate. 

Third, the present result does not exclude that the key element for diffusion anomaly is the 
positive correlation of density and energy, instead of entropy. While in water and the present 
model the LDL phase has lower energy and entropy than the HDL, in the isotropic potentials with 
water-like anomalies the LDL phase has lower energy but higher entropy than the HDL. These 
considerations support the idea that the mechanism of anomalies in isotropic potentials is different 
from that of water. 9 

Interestingly, here we also observe that the lines of constant Dii resemble the melting line of 
bulk water. At low temperatures, we find the locus of density maxima, which marks another well- 
known water anomaly. We discuss how this locus is related to the locus of expansivity minima. 

At lower T, we find subdiffusive behavior Ar 2 ~ t a , as seen in experiments of hydration wa- 
ter— and simulations of confined water.— iI2Iril£ Our results are well described by a = 0.7, be- 
tween 0.96^ and 0.45 ± 0.05-^ of previous calculations. Previous works proposed that subdiffu- 
sion is a consequence of the heterogeneities in water-interface interaction. Here this rationale does 
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not apply and we relate the subdiffusion to the increase of correlation among water molecules at 
low T due to the full development of the HB network. 

By further decreasing T , we find density minima, as seen in experiments 69 and MD simula- 
tions.— 1 ^- These minima occur within the subdiffusive part of the phase diagram, therefore where 
translational motion is strongly hampered and glassy behavior is incipient. Nevertheless, the HB 
network within this region is still dynamically evolving, with increasing correlation time. 16 In 
particular, the HB correlation time is about 4 ns for the majority of the subdiffusive region and 
increases, exceeding our simulation times of the order of 15 ms, only at about P < 0.02 GPa and 
T < 170 K. 

Previous works related the density minima to the saturation of a network of molecules with a 
random tetrahedral local arrangement. 71 However, this rationale apply to our model only for the 
case in which the cooperative component of the hydrogen bonds is zero. When the cooperative 
component is larger than zero, as expected in real water, 25 our calculations show that the minima 
are due to a reduction of the number of HBs, as a consequence of the reorientation of molecules 
for a better matching of local order. The high energy cost of this rearrangement is the cause of the 
large slowing down of the HB dynamics near the state points where the density minima occur. 
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